Short-range type critical behavior in spite of long-range interactions: 
the phase transition of a Coulomb system on a lattice 
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One- to three-dimensional hypercubic lattices half-filled with localized particles interacting via the 
long-range Coulomb potential are investigated numerically. The temperature dependences of specific 
heat, mean staggered occupation, and of a generalized susceptibility indicate order-disorder phase 
transitions in two- and three-dimensional systems. The critical properties, clarified by finite-size 
scaling analysis, are consistent with those of the Ising model with short-range interaction. 

PACS numbers: 64.60.Fr, 05.70.Jk, 71.10.-w, 02.70.Uu 



The problem under which conditions Coulomb glasses 
exhibit genuine phase transitions has been under con- 
troversial debate for two decades P, S S 0|- The 
significance of the kind of static disorder is unclear yet 
Is, 4J. In this context, the relation to the Ising model 
with short-range interaction is of great interest. It would 
be very useful to know how replacing its nearest-neighbor 
coupling by an "antiferromagnetic" long-range Coulomb 
interaction modifies the critical behavior of a system 
without static disorder. Due to the interplay of long- 
range interaction and frustration, one might alterna- 
tively expect this model to belong to different univer- 
sality classes: In case of a sufficiently slowly decaying 
ferromagnetic power-law interaction, including the l/r- 
interaction, mean field behavior is obtained But the 
condensation in a three-dimensional model of an elec- 
trolyte exhibits Ising-like critical behavior, possibly due 
to efficient screening 0- 

To decide the question, we numerically study lattices 
half-filled with localized particles interacting via the long- 
range Coulomb potential. 



H = 



(n,; - 1/2) - 1/2) 



(1) 



Here ria G {0, 1} denote the occupation numbers of states 
localized at sites within a d-dimensional hypercube of 
size L'^. Elementary charge, lattice spacing, dielectric 
and Boltzmann constants are all taken to be 1. Neutral- 
ity is achieved by background charges -1/2 at each site. 

For reducing finite-size effects, we impose periodic 
boundary conditions for o? = 1 and 2, using the mini- 
mum image convention (^. For d = 3, the same approach 
would give rise to an unphysical feature: The ground- 
state would be a layered arrangement of charges instead 
of the expected NaCl structure in case L is a multiple 
of 4 Hence, for d = 3, we consider the sample to be 
surrounded by eight equally occupied cubes. 
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Our numerical investigations have been performed by 
means of algorithms based on the Metropolis procedure 
IE IE ^1 ■ Such simulations are very expensive not only 
for the long-range interaction, but also for correlation 
times diverging close to phase transitions and as temper- 
ature T vanishes. Therefore it is necessary to adapt the 
dynamics to the situation that the simultaneous change 
of rii for certain clusters is necessary for overcoming 
passes in the energy landscape. To our knowledge, a 
procedure similar to the Swendsen-Wang algorithm |ll| 
is not available for the interaction type considered here. 
Hence we modified the dynamics "by hand" to take into 
account several kinds of processes: one-electron exchange 
with the surroundings, one-electron hops over (T depen- 
dent) restricted distance, and two-electron hops simul- 
taneously changing the occupation of four neighboring 
sites. 

At high T, we use the original Metropolis method 
But at low T, we take advantage of a hybrid procedure 
much accelerating the computations. It connects the di- 
rect evaluation of weighted sums over states within a low- 
energy subset of the configuration space with Metropolis 
sampling of the complementary high-energy subset 10]. 

We now turn to the qualitative behavior of specific 
heat, order parameter, and susceptibility. The specific 
heat c was obtained from energy fluctuations utilizing 
c =: {{H^) - (iJ)2)/(T2 L"^) . The T and L dependences 
of c are presented in Fig. 1 for dimensions d = 1 to 3. 
This graph shows that, for d = 2 and 3, sharp peaks of in- 
creasing height evolve within a small T region as L grows. 
Away from the peaks, within the T intervals presented, 
c is almost independent of L. But for d = 1, there are 
only broad, rounded peaks with L independent height - 
a logarithmic T scale is used for d = 1, in contrast to the 
linear scales for d — 2 and 3 which display far smaller T 
intervals. For d = 1, finite size effects are restricted to 
low T where the reliability bound decreases with L. 

Hence, according to the behavior of c(T, L) a phase 
transition likely occurs for d = 2 and 3, in agreement 
with results of lattice gas simulations for d = 3 0, 0| ■ 
However, for d = 1, in spite of the long-range interaction. 
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FIG. 1: Temperature dependences of the specific heat, c{T), 
for dimensions d = 1 to 3 as obtained from simulations of 
samples of L'' sites, d = 1: L = 100 (+), 280 (o), and 700 
(x); d = 2: L = 20 (+), 34 (o), and 58 (x); d = 3: L = 8 (+), 
12 (o), and 18 (x). Only a part of the data points forming the 
curves is marked by symbols. The error bars are considerably 
smaller than the symbol size. 



there seems to be no phase transition at finite T. 

Analogously to an antiferromagnet, the order inher- 
ent in a charge arrangement rii can be characterized by 
means of the staggered occupation ai relating to a NaCl 
structure. For example, if d = 3, it is given by 



= (2n, - 1) . (-iy^+y^+^' 



(2) 



where Xi, i/i, and Zi denote the (integer) components of 
r;. Thus we measure the mean absolute value of the 
staggered occupation {\a\) as order parameter. 

T and L dependences of (|(t|) are shown in Fig. 2. For 
d = 1, a rapid decrease of {\a\) with increasing T oc- 
curs already clearly below the temperature of maximum 
c (same T scales in Figs. 1 and 2). This marked decrease 
shifts to lower T with increasing L. For = 2, a qualita- 
tively different behavior is found: decreases rapidly 
just in that T region where the peak of c(T) evolves. 
The T interval of rapidly diminishing {\a\) shrinks as L 
rises, another indication of the phase transition. This 
interpretation is confirmed by an extrapolation i — > oo: 
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FIG. 2: Temperature dependence of the mean absolute value 
of the staggered occupation, {|cr|)(T), for d = 1 and 2. For 
the meaning of the symbols o, and x see caption of Fig. 
1; • marks the extrapolation L — > oo explained in the text. 
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FIG. 3: Temperature dependence of the susceptibility, xiT), 
related to the staggered occupation a by Eq. (3), for d = 1 
and 3. For the meaning of the symbols see caption of Fig. 
1. Additionally, for d = 1, the inset shows how Tmax, the 
temperature of maximum x{T), depends on sample size L. 



Assuming {\a\){T,L) = {\a\) {T, oo) + A{T) / L'^''^ with A 
independent of L, an almost sharp transition is obtained 
from data for L — and 58, although this extrapola- 
tion is of limited accuracy in the immediate vicinity of 
the transition. For d = 3, the results (not shown here) 
qualitatively resemble our findings for d = 2. 
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FIG. 4: Temperature dependence of qd, related to the Binder 
cumulant, in the close vicinity of the transition for d = 2 and 
3. For definitions see text. With increasing modulus of the 
slope, the curves refer to L = 16, 24, 34, 48, 68, and 96 for 
d = 2, and to L = 8, 10, 12, 14, 16, and 18 for d = 3. For 
clarity, data points in the intersection regions are omitted. 

The generalized susceptibility x related to the order 
parameter {\a\) is given by 

X^L''{{a')^{\a\r)/T . (3) 

Fig. 3 shows T and L dependences of x- It confirms the 
conclusions drawn from (|(t|)(T, L): On the one hand, for 
d = 1, a broad peak of x(^) evolves with increasing L 
where Tmaxj the temperature of maximum x, decreases. 
As the inset demonstrates, T^ax can be approximated 
by a/\n{bL) with constants a and b. Hence, Tmax hkely 
vanishes as i ^ cx) so that there seems to be no phase 
transition for d = 1 at finite T. On the other hand, for 
d = 3, as L rises, a narrow peak grows in just that T 
region where c(T, L) has such a feature, a further indica- 
tion of the phase transition. For d — 2, x{T, L) behaves 
qualitatively similar to the results for d = 3. 

The quantitative evaluation of our simulation data con- 
sists in a finite-size scaling analysis 0, For this 
aim, we first consider q2 = — ln(l — ((T^}^/(cr'*)) and 
93 = -tan(7r(l - l.h {a^ I {f^^))) for d = 2 and 3, re- 
spectively. These quantities are directly derived from the 
Binder cumulant 1 — (cr^)/(3(CT^)^) and exhibit similar 
features. The qd{T) have small curvature in the vicinity 
of the transition what alleviates a precise interpolation. 

Fig. 4 shows qd{T). For d — 2, there clearly is a com- 
mon intersection point of the curves for different L at the 
critical temperature Tc.2- But for d = 3, only a tendency 
towards such a behavior is seen. However, the system- 
atic corrections to scaling can to a large extent be taken 
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FIG. 5: Relation between asiL) and Tc,3(L) for 53,0 = -0.174. 
With increasing as, the points refer to L = 8, 10, 12, 14, 16, 
and 18. The error bar represents our extrapolation L 00. 

into account in a simple way: We define a size-dependent 
critical temperature by qd{Tc^diL), L) = g^ o where the 
qdfi are appropriate constants, fixed below. Scaling of 
qd{T) curves for different L with respect to T — Tca{L) 
yields very good data collapse. (Referring to t = in- 
stead of T = Tc^d{oo) considerably reduces the infiuence 
of deviations from scaling on the values for the critical 
exponents Thus we assume qd to depend only on 

t = ad{L){T — Tc^diL)) in the vicinity of the transition. 

We approximate qdit) by polynomials of third degree, 
qd,o + t + bdt^ + Cdt^- By regression studies we first ad- 
justed the L independent parameters bd and Cd of the 
ansatz, and then determined the ad{L) values. The lat- 
ter were analyzed by means of power law fits including 
various L intervals, where according to finite-size scaling 
ad{L) oc L^^'^ was presumed. For consistency, the mean- 
square deviation of these fits must be understandable as 
resulting from random errors alone. Tab. I presents the 
most precise results for u, which were obtained from the 
fits safely fulfilling this requirement. 

In obtaining Tc,d{L) from qd{TcMiL),L) = qd,o, a de- 
viation S of qdfi from the L ^ 00 limit of qd{T,2L) = 
qd(T, L) gives rise to a contribution (x 5/ad{L) to Tc^d{L). 
We chose qd,o so that this term vanishes: 92,0 = 1.8933 
and (73 = —0.174. The remaining higher order correc- 
tions in Tc^d{L) originate from imperfection of finite-size 
scaling. Comparing several empirical approximations, 
we observed that, over a wide L range, they are al- 
most proportional to ad{L)~'^, see Fig. 5. Correspond- 
ing extrapolations yield the following values of Ted (00): 
0.10308 ± 0.00002 and 0.12881 ± 0.00003 for d = 2 and 
3, respectively. The confidence intervals include the 3<t- 
random errors and cautious estimates for the systematic 
uncertainty of the extrapolation L — > 00, see Fig. 5. 

The analysis of c(T,L), (|ct|)(T,L), and x(T, L) was 
performed similarly to the evaluation of qd{T,L): We 
considered Inc, ln(|cr|), and Inx as functions of t and L. 
For not too large as L ^ 00, scaling implies that each 
of these quantities is decomposable into a sum of two 
functions depending only on t and L, respectively. How- 
ever, for the L regions considered here, this hypothesis 
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FIG. 6: Size dependence of the value of (|cr|) at T^^diL) defined 
in the text. Error bars are considerably smaller than symbol 
size. The dashed lines represent the fits given in Tab. I. 



Quantity d L region Coulomb s.-r. Ising 



a/v 
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24 


- 96 


-0.037(44) 


(In) 


/3/u 


2 


48 


- 96 


0.131(6) 


1/8 


■y/u 
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48 


- 96 


1.735(24) 
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34 


- 96 


1.021(33) 
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ajv 


3 


6 - 


- 18 


0.09(6) 


0.1761 [25] 


P/v 


3 


12 


- 18 


0.499(16) 


0.5181[4] 


-y/u 


3 


14 


- 18 


1.974(23) 


1.9638[8] 


h> 


3 


10 


- 18 


0.635(10) 


0.6297[5] 



TABLE L Finite-size scaling results for the critical exponents 
of specific heat, staggered occupation, susceptibility, and cor- 
relation length, a, 13, 7, and v, respectively. To retain nu- 
merical precision, we mostly present exponent ratios instead 
of the exponents themselves. Values for the Ising model with 
short-range interaction fll . [Tst are included for comparison. 
Parentheses and brackets give 3cr-random errors and total er- 
rors, respectively, referring to the last given digit of the value. 



took into account a background contribution presuming 
c{Tc{L),L) = f{L) = a + bL"'^". However, we treated 
f{Li) and /(i2) at fixed finite Li and L2 as adjustable 
parameters instead of a and b to avoid numerical prob- 
lems for almost logarithmic behavior of c(L) (small \a\). 
Results for a/iy obtained this way are included in Tab. I. 

Our values in Tab. I have to be regarded as effective 
exponents. Due to the finiteness of L, tiny systematic er- 
rors are certainly present, presumably the more relevant 
the smaller the exponent value. Unfortunately, our data 
set is not sufficient for a convincing estimate of them. 
However, even if only random errors are considered, our 
values comply with the Widom relation, 2 = a + 2/3 + ^, 
and the hyperscaling relation, 2 — a = dv. 

For comparison. Tab. I includes also values for the Ising 
model with nearest- neighbor interaction 0, It is 

surprising that, in spite of the differing characters of the 
interactions, the critical exponent data obtained here are 
very close to these values: In particular, for 7/1/ and v. 
the agreement is perfect within numerical accuracy. The 
slight deviations in a/iy and P/iy for d — 3 presumably 
arise from a too small sample size. Note: obtaining these 
values indirectly, via Widom and hyperscaling relation 
from 7/ J/ and v, yields again perfect agreement. 

Concluding, in spite of the long-range interaction, the 
Coulomb system described by Eq. (1) seems to belong 
to the same universality class as the Ising model with 
short-range interaction. This suggests that the lattice 
Coulomb-glass model might have the same critical prop- 
erties as the random-field short-range Ising model. 

We thank H. Eschrig, M.E. Fisher, B. Kramer, T. Nat- 
termann, M. Richter, M. Schreiber, and T. Vojtafor help- 
ful discussions and literature hints. 



proved to be well fulfilled only for In^. In the cases of 
Inc and Incr, there is a clear tendency towards such a be- 
havior, but small deviations cannot be neglected. Thus 
we approximated Inc, ln(|cr|), and Inx by polynomials in 
t of third order, taking advantage of universalities in the 
coefficients as far as possible. This regression provides 
precise values for the observables a.t t — 0. Simulta- 
neously, we obtained the confidence intervals taking into 
account the uncertainties in the individual measurements 
of the observables and in the Tc^d{L) values. 

The interpolation results for {\a\) and x were ana- 
lyzed by means of power law fits, where proportionality 
to L'^/" and L'^/^ , respectively, was presumed. These 
studies were performed analogously to the determination 
of V. Tab. I presents the values of the exponent ratios 
(3/u and 'y/v. The high quality of these power law fits is 
illustrated by Fig. 6 presenting (|(t|)(Tc(L), L). 

Compared to the study of (|ct|) and X) the analysis of 
c is more difficult: Exponent values obtained from power 
law fits converge only slowly with increasing L, and the 
mean-square deviations remain too large. Therefore we 



[1] 

[2] 
[3] 

[6: 

[7 

[9: 

[10 

[11 

[12 
[13 



26, 2883 (1993). 
Rev. Lett. 71, 3335 



J.H. Davies, P.A. Lee, and T.M. Rice, Phys. Rev. Lett 
49, 758 (1982). 

T. Vojta, J. Phys. A: Math. Gen. 
E.R. Grannan, C.C. Yu, Phys. 
(1993). 

T. Vojta, M. Schreiber, Phys. Rev. Lett. 73, 2933 (1994). 
A. Dfaz-Sanchez, M. Ortuno, A. Perez-Garrido, 
E. Cuevas, phys. stat. sol. (b) 218, 11 (2000). 
E. Luijten and W.J. Blote, Phys. Rev. Lett. 89, 025703-1 
(2002), and refs. therein. 

E. Luijten, M.E. Fisher, and A.Z. Panagiotopoulos, Phys. 
Rev. Lett. 88, 184701-1 (2002). 

N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, 
A.M. Teller, E. Teller, J. Chem. Phys. 21, 1087 (1953). 
A. Mobius, P. Thomas, J. Talamantes, and C.J. Adkins, 
Phil. Mag. B 81, 1105 (2001). 

A. Mobius, P. Thomas, Phys. Rev. B 55 , 7460 (1997). 
R.H. Swendsen and J.S. Wang, Phys. Rev. Lett. 58, 86 
(1987). 

R. Dickman and G. Stell, |cond-mat /9906364| 

K. Binder, E. Luijten, M. Miiller, N.B. Wilding, and 

H.W.J. Blote, Physica A 281, 112 (2000), and refs. 



5 

therein. (Wesley, Reading, 1992). 

[14] N. Goldenfeld, Lectures on Phase Transitions and the [15] M. Hasenbusch, Int. J. Mod. Phys. C 12, 911 (2001). 
Renormalization Group, Frontiers in Physics, Vol. 85, 



